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SUMMARY 


This report contains our studies on applications of the finite-element 
approach to transonic flow calculations, and it includes comparisons between 
different discretization techniques of the differential equations and boundary 
conditions. Transonic flow calculations can be divided into two main categories: 
type-sensitive methods and type- insensitive methods. Finite-element analogs of 
Murman’s mixed-type finite-difference operators for small-disturbance formula- 
tions are constructed, with different strategies used in the subsonic and super- 
sonic regions. In the supersonic region, no upstream effect is allowed; 
blending elements are introduced between different regions. On the other hand, 
as an example of type-insensitive methods, the time-dependent (unsteady) approach 
(using finite differences in time, finite elements in space) is examined. The 
elliptic methods, where the transonic equation is cast into Poisson’s form with 
the nonlinear terms as a driving force, provide another example. The report is 
concluded with a general shock-fitting procedure based on discontinuous shape 
functions and with possible extensions to full potential equations. 


INTRODUCTION 


Computations of steady transonic flow can be formulated in terms of 
either the Euler equations or the velocity potential equation; but, regard- 
less of which formulation is chosen, these computations generally rely on one 
of two basic iterative procedures. The first procedure involves integrating 
a set of hyperbolic equations (in time) until a steady state is reached, while 
the second approach makes use of relaxation techniques. 

The hyperbolic procedure is frequently formulated in terms of the com- 
plete unsteady Euler equations, but ether hyperbolic forms of the equations 
of motion have also been used. The hyperbolic procedure is attractive 
because a converged solution, which includes both subsonic and supersonic 
regions, can be obtained without making any explicit consideration of the 
mixed character of the flow field. Unfortunately, the convergence to the 
steady state has proved to be quite slow. 

By contrast, the relaxation procedures cannot produce converged solutions 
unless special local (spatial) discretization procedures that account for 
the mixed elliptic-hyperbolic nature of the flow field are used. When the 
local character of the flow field is properly accounted for, however, the re- 
laxation procedures converge much faster than the hyperbolic ones. The choice 
of the most appropriate type of spatial discretization must take into account 
their differences. 

Solutions of the Euler equations in the transonic regime have been ob- 
tained by means of standard finite-difference procedures (based on the 
Lax-Wendroff or the McCormack schemes) as well as finite-volume techniques. 
Solutions cf the velocity potential equation have generally relied on the use 
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of type-dependent, finite- difference schemes cf the kind first introduced to 
numerical transonic flow calculations by Murman and Cole (ref. 1). Recently, 
accelerated iteration procedures for finite-difference calculations using 
fast Poisson solvers (elliptic methods) proved to be very efficient. Martin 
and Lomax (ref. 2) used elliptic methods, in the form of a system of first- 
order equations, for cases of small disturbance, while Jameson (ref. 3) solved 
the full potential equation. 

In this report, we examine the application of the finite- element approach 
to transonic flow calculations. We consider hyperbolic, mixed- type, and 
elliptic methods. The appeal of finite-element procedures is twofold. First, 
finite-element procedures are capable of accurately and efficiently enforcing 
boundary conditions, even when the boundaries are geometrically complex. 

(The application of boundary conditions in finite-difference schemes becomes 
very difficult when the boundaries are complex in shape.) Second, finite- 
element procedures reduce the number of grid points (or elements) required 
to achieve a solution of a desired accuracy through the use of efficient, 
higher-order shape functions or mixed finite- element methods. We note that, 
although both of these advantages are important in two-dimensional flows 
(with which we are concerned in this report), they become crucial in three- 
dimensional calculations. 

We have carried over from finite-difference methods as much understanding 
of numerical transonic techniques as possible. More specifically, we have 
made the basic assumption that techniques that are successful in transonic flow 
regimes when using finite-difference methods should also be successful when 
using finite-element techniques, and likewise for techniques that fail. As a 
result, our primary emphasis has been on solving the small-disturbance equation. 
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not the full velocity potential equation. Furthermore, we have been satisfied 
with enforcing boundary conditions on a straight line (in accordance with 
small-disturbance theory) . The use of finite-element procedures on this sim- 
plified problem as an end in itself is not justified; finite-difference methods 
are undoubtedly better suited to the small-disturbance problem. Nevertheless, 
we believe that the feasibility of developing finite- element procedures that 
are capable of handling mixed, elliptic-hyperbolic flow fields can best be 
demonstrated in this simpler environment. Extension to the more complicated 
problems for which finite-element methods are better suited should be simpli- 
fied after the small- disturbance problem has been completed. 

Much of our present effort, therefore, has been spent on solving the 
small-disturbance equation in a rectangular domain with simplified boundary 
conditions, by means of rectangular elements and linear shape functions. This 
simplification not only makes the finite-element procedure easier to apply, 
but also brings it parallel with finite-difference procedures, which have 
proven to be successful in this simplified problem. 

From this "jumping-off" point, we proceed to higher-order shape functions 
with a parallel review of higher-order-accurate, finite-difference procedures 
so that effects that are related to improved accuracy can be separated from 
those that are derived from the basic differences between finite-element and 
finite-difference techniques. Finally, using elements with curved boundaries, 
we consider the solution of the full potential equation and its extension to 
nonrectangular domains. 

One of the most prominent differences between finite differences and 
finite elements that persists even in the simplified problem described above 
is that the matrices associated with finite-element schemes are generally 
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solved by direct methods rather than by the relaxation techniques that are used 
in finite-difference methods. Economically, iterative techniques are usually 
more suitable for large, sparse matrices, whereas direct methods are more 
suitable for more moderately sized matrices. As noted above, _J:he larger ele- 
ments that can be used for finite elements result in matrices more moderate in 
size, so direct inversion is favored. Later in this report we consider whether 
it is more economical to relax or to invert directly a given matrix system, but 
we note here that there is mathematical basis for selecting one scheme over the 
other. For example, because of the unique properties that are transmitted to 
the matrix when the equation changes type (from elliptic to hyperbolic) , re- 
laxation methods could fail while direct methods would succeed. That is, the 
relaxation procedure might diverge even in cases where the inverse of the 
matrix actually exists. 

One of the most important features of transonic flows is shock waves. 

Like finite-difference methods, shocks are either captured or fitted. In this 
report we discuss a general shock-fitting procedure for finite-element calcula- 
tions and we use a simple version in our computations. 

TRANSONIC SMALL-DISTURBANCE THEORY 

The assumption of the existence of a velocity potential, along with re- 
striction to small disturbances, greatly simplifies the transonic flow problem 
while, at the same time, it retains all the fundamental nonlinear, mixed- type 
mathematical properties that are characteristic of transonic flows. Despite 
its relative simplicity, the small-disturbance equation is capable of ade- 
quately describing the transonic flow field around many configurations of 
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practical Interest, and it has been used in numerous engineering applications. 
In its small-disturbance form, the transonic velocity potential equation is 

(K- (j) )(f) + (J) = 0 , (1) 

^xx ^yy 

where ^ represents the perturbation velocity potential, K represents a 
transonic similarity parameter, and x and y refer to a Cartesian'^oordi- 
nate system (see fig. 1). The boundary condition to be applied in the vicinity 
of the body is 

(f)y = f^(x) - a , (2) 

where represents the body shape and a signifies the angle of attack. 

Consistently with the small-disturbance approximation, this boundary condition 
may be applied on the line y = 0 . The airfoil surface boundary condition is 
generally complemented by an analytic expression for the far field, which re- 
presents the combined effects of a doublet and a vortex. The formulation of 
the problem is summarized in figure 1. 

In the small-disturbance equation, the sonic line is given by the point 
where = K , and, by inspection of equation (1) , we see that when 

d) = K (sonic condition), d) vanishes. This latter condition is enforced 
^x yy 

explicitly in many numerical schemes. Once (|)^ becomes greater than K , 
the equation changes type and becomes hyperbolic. In the hyperbolic (super- 
sonic) region, the characteristic directions are 

(dx/dy)^ = - K , (3) 
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while the shock relations for the small-disturbance velocity potential for- 
mulation are 


<K - (j)^> = - (dx/dy)^ 


( 4 ) 


and 



(5) 


where < > and ([ |] signify the average and the jump across the dis- 

continuity, respectively. 

Related Work on the Small— Disturbance Equation 

In conjunction with our discussion of the small-disturbance equation, we 
mention some related work. Some of this work is related to information that 
can be used to understand the iterative procedure, some is concerned with de- 
veloping more economical iterative procedures, and some describes appropriate 
limitations on the small-disturbance approximations that support the transonic 
small-disturbance equation. 

One interesting study is the work of Sichel (ref. 4), who considered the 
effects of viscosity on transonic flows. In his work, Sichel used the viscous 
transonic equation: 


{}) (j) - (j) = V(f) 

X XX yy XXX 


( 6 ) 
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Although his work points out some important limitations of inviscid transonic 
small-disturbance theory, we wish to reference here the direct parallel be- 
tween the physical viscosity term in the viscous transonic equation and the 
artificial (or numerical) viscosity that is used in numerical transonic flow 
studies. We discuss artificial viscosity in later sections. 

Another work of interest is Landahl's investigation (ref. 5) of the un- 
steady transonic equation: 

(K - (}) + ac{)^)<|) + (t> = *. + • - (7) 

^t XX ^yy ^xt ^tt 

This equation can be construed as representing one of the complete unsteady 
solution procedures described in the Introduction. A second interesting aspect 
of this equation, however, is the so-called low-frequency form of this equation, 
which is obtained by neglecting the high-frequency terms and (J)^ as 

follows: 


(K - (|) )(j) + (f) = ^ • (8) 

^x XX yy ^xt 

Equation (8) is closely related to the relaxation procedure used to solve the 
transonic small-disturbance equation. 

Finally, we should enumerate some of the limitations of the transonic 
small-disturbance equation. The approximations upon which the transonic small- 
disturbance equation are based generally break down near the leading edge of a 
transonic wing. Such leading edges are, for engineering reasons, generally 
blunt, so the flow must turn as much as 90° from its original direction. Turns 
of this magnitude are not allowed in small-disturbance theory. Keyfitz, Melnik, 
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and Grossman (ref. 6) have given more complete consideration of the problem of 
the blunt leading edge. Fortunately, ignoring the details of the leading edge 
region allows acceptable engineering accuracy to be achieved in the remainder 
of the flow field. 

Sirovich and Huo (ref. 7) have tested the validity of the transonic small- 
disturbance equation in the vicinity of the sonic line, while Landau (ref. 8) 
and Guderley (ref. 9) have discussed the details of the flow in the intersection 
between the sonic line and the shock wave. The intersection of a normal shock 
wave with a curved surface has been discussed by Zierep and Oswatitsch (ref. 10), 
who determined the character of the solution near this singularity. 

One transonic phenomenon that we do not consider in this report, but must 
at least mention, is the effect of viscosity in transonic flow regimes. Vis- 
cosity is important in shock-wave boundary- layer interactions and in the 
trailing edge region. Either of these regions can generate local-separation 
bubbles, which substantially alter the flow from its unseparated, inviscid 
state. Some specific works that discuss methods for including these viscous 
effects include those by Melnik and Grossman (ref. 11). 

FINITE DIFFERENCES 
Unsteady Approach 

Magnus and Yoshihara (ref. 12) obtained numerical solutions of the Euler 
equation in the transonic region, using a Lax-Wendroff finite-difference scheme 
(with artificial viscosity) marching in time to the steady state. For small 
disturbances, Magnus and Yoshihara used the following hyperbolic system of 
equations : 
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( 9 ) 


U = (K - U)U + V , 
t. X y 


V 


t 


= -V + u . 
X y 


( 10 ) 


The calculations were too lengthy and expensive, and, hence, this method 
was abandoned. 


Murman's Fully Conservative Scheme 

One of the most popular techniques for solving the transonic small- 
disturbance equation (TSDE) is Murman^s fully conservative, type- dependent , 
finite-difference scheme, or variants of it (see ref. 13). As indicated in 
figure 2, Murman’s scheme is characterized by four distinct operators: an 

elliptic operator E for subsonic regions; a hyperbolic operator H for 
supersonic regions; a parabolic operator P for points on (or near) the 
sonic line; and a shock-point operator SPO for enforcing the jump conditions 
across the shock. The elliptic operator is based on a second-order-accurate, 
central-difference formula, while the hyperbolic operator is obtained from a 
first-order-accurate, backward-difference representation. The remaining two 
operators, the parabolic operator and the shock-point operator, represent 
blending elements for grid points on the boundaries between elliptic and hy- 
perbolic regions. 

The four finite— difference operators described above can be used to con- 
vert the continuous, partial differential equation (1) into a discrete system 
of algebraic equations which describe the behavior of the solution at a 
fixed set of points in the flow field. The system of equations generated 
by this discretization is generally solved by a line-relaxation algorithm 
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which iteratively sweeps the matrix, line-by-line, from left to right. 

During the relaxation process, the solution is over-relaxed (o) > 1) in the 
elliptic region, but in the remainder of the field it is usually mildly under- 
relaxed (w < 1) , Since this discretization scheme is only first-order 
accurate in the supersonic region, it is first-order accurate overall. Be- 
cause of the relative simplicity of the scheme, the low accuracy can be 
offset by using a finely divided mesh. 

We note that, in the far field, the doublet strength due to the flow 
around the body is given by 

D = // dx dy + . (11) 

As can be seen, the doublet strength depends on the solution itself. The 
double integral is generally updated at selected intervals during the 
iterative procedure so that when the solution converges, the far-field 
representation is intimately tied to the numerical solution (and conversely) . 
The use of this analytical representation of the solution far from the 
airfoil greatly decreases the domain included in the computation. 

A series of sample calculations, which have been obtained from Murman's 
fully conservative scheme and which indicate the influence of various factors 
in the scheme on the final solution, are given in figures 3A through 3C. 


We also note that Cheng and Hafez (ref. 14) have shown that the far-field 
behavior near the boundary may be fitted by using a least-squares technique; 
thus, performing the double integration over the entire flow field is no 
longer necessary. 
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The base case, which is designated as "Case Al" in the figure, includes all 
the factors in the Murman scheme, as described above. Unless otherwise noted, 
these calculations are for a uniformly-spaced grid system in both the x 
and y directions, with 10 grid points on the airfoil surface, 10 points up- 
stream of the leading edge, and 10 points downstream of the trailing edge. 

The first comparison (fig. 3A) demonstrates the effect of the shock- 
point operator on the final results. As shown, the shock-point operator 
allows a much more rapid velocity change across the shock (compare Cases Al 
and B1 in fig. 3B). The comparison between Cases Al and Cl shows the effect 
of incorporating the far-field solution (Case Al) instead of enforcing a 
homogeneous boundary condition (4> = 0) at the same point (Case Cl) . Case D1 
shows the effect of placing the leading and trailing edge points halfway be- 
tween grid points, compared to the effect of placing grid points on the leading 
and trailing edge points. The results in figure 3C show that the apparent 
shape of the airfoil is altered when the positions of the leading and trailing 
edge points are shifted with respect to the grid system. This effect is 
amplified in the present case by our use of a relatively coarse grid system. 
Also in figure 3C, we show the effect of halving the grid size. Cases D1 and 
HI show the higher degree of accuracy that can be obtained with the finer 
grid system. 

Among other considerations that are directly related to numerical methods 
are those techniques concerned with accelerating Murman' s relaxation solution 
so that convergence can be obtained more rapidly. One such approach is the 
use of extrapolation techniques, such as those reported by Cheng and Hafez 
(ref. 15) and Caughy and Jameson (ref. 16). Extrapolation techniques attempt 
to accelerate the painfully slow iterative procedure by obtaining numerical 
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estimates of the dominant eigenvalues in the finite-difference matrix and by 
using these eigenvalues to extrapolate the Iteration to a level much closer 
to convergence. Such extrapolation techniques have exhibited time savings of a 
factor of 4 to 5 over the more conventional line-relaxation procedures. 

A second acceleration device, which has proven very effective for 
finite-difference calculations, is one that replaces the nonlinear transonic 
equation with a Poisson equation having the nonlinear term as a right-hand-side 
term. The matrix corresponding to the constant-coefficient elliptic (left- 
hand-side) operator is then solved by a direct inversion technique. Based 
on this interim solution, the nonlinear (right-hand-side) term is updated, and 
the "fast-solver" is again employed. With this technique, convergence has been 
very rapid (generally less than 10 of these major iterations), provided the 
supersonic region is treated properly. If no special treatment is applied 
for the supersonic region, the procedure fails. Very impressive results (in 
terms of the amount of computer time required) have been reported by both 
Martin and Lomax and by Jameson (see refs. 2 and 3). In their studies, 
slightly different techniques were used for treating the supersonic region. 

These are discussed more fully in the following section, A Fully Conservative, 
Second-Order Scheme for Finite Differences. 

Two additional studies that are concerned with developing more rapid com- 
putational procedures are also in progress. Ballhaus and Steger (ref. 17) 
and Jameson (ref. 3) are using alternating-direction implicit methods, based on 
an efficient matrix factorization technique for solving the Euler equations. 

The second approach is a multi-grid technique being developed by Brandt and South 
(ref. 18). This technique solves the transonic equation on a series of grid 


13 



spacings which vary from fine to coarse and fine again. The basis of this pro- 
cedure is to increase the speed of the iterative relaxation process by 
diminishing the error in a different segment of the frequency spectrum of the 
matrix with each distinct grid system. 

Another acceleration technique is the use of discrete, shock-fitting pro- 
cedures, which have been reported by Cheng and Hafez and by Yu and Seebass 
(see refs. 15 and 19). Shock-fitting really has no capability for increasing 
the speed of the relaxation process; however, it indirectly achieves such £.n 
improvement. Shock-fitting is really concerned with improving the spatial 
discretization so that shock waves can be handled with a relatively coarse mesh 
(compared to the extremely fine grid that must be used when the shock is not 
fitted). The coarse grid that is used with shock-fitting decreases the amount 
of computational time by reducing the number of unknowns (grid points) to be 
computed. 

A Fully Conservative, Second-Order Scheme for Finite Differences 

As Indicated in the previous section, Murman's fully conservative finite- 
difference scheme is accurate only to the first order. The low accuracy occurs 
because upstream effects are not allowed in the supersonic region, and, hence, 
one-sided differences must be used in that region. The numerical scheme is 
therefore (formally) less accurate in supersonic regions than in subsonic 
regions, where central differencing is employed. To make the accuracy in the 
supersonic region comparable with that in the subsonic region (thus making 
the entire calculation second-order accurate), we may use either of two methods. 
We can include more grid points in the computations for supersonic points, or we 
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can turn to a Hermitian scheme in which both (p and its derivatives are stored 
and used at each grid point. 

One reason for considering a finite— element scheme is the premise that, 
through the use of a more accurate local representation of the solution, the 
total number of grid points can be decreased, and, thus, the computer storage 
requirements, as well as the central processor time requirements, can be 
reduced. Such goals are particularly urgent for three-dimensional transonic 
flows. With the goal of improved accuracy in mind, we consider first the 
requirements that are imposed on a finite-difference scheme when it is extended 
to higher accuracy. It seems reasonable to expect that those problems encoun- 
tered in a second-order-accurate, finite-difference scheme would also be present 
in a finite-element scheme, so their solution should give some indication of 
whether a given finite-element technique will be successful. 

Murman and Cole suggest an implicit, second-order-accurate, backward 
difference scheme as an alternate to the first-order-accurate, backward dif- 
ference scheme that was used in the supersonic flow region. They hinted in 
their paper (see ref. 1), and it has been verified since, that calculations 
employing this second-order-accurate, backward difference scheme sometimes 
diverge during the relaxation iteration. Because of this difficulty, nearly 
all velocity potential calculations have been made with schemes accurate only 
to the first-order. 

More recently. Warming and Beam (see ref. 20; see also Martin’s work, 
ref. 21) have extended Murman' s fully conservative scheme to the Euler 
equations. Although they were primarily concerned with the Euler equations. 
Warming and Beam suggested (but did not test) a second-order-accurate, back- 
ward difference scheme for the transonic small-disturbance equation, including 
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a shock-point operator and a parabolic point. We note that Murman’s, as well 
as Wanning and Beam’s, second-order backward difference scheme is both dis- 
sipative and dispersive compared to Richtmeyer’s standard (non-dissipative) 
scheme (see figs. 4A and 4B) . The important questions concerning these 
second-order schemes, especially in our finite-element context, are why 
is the second-order-accurate scheme of Murman so much less reliable than 
his first-order scheme and will the alternative proposal of Warming and Beam 
alleviate this problem? 

We begin by considering the latter of these two questions. We incorporated 
the second-order-accurate scheme of Warming and Beam into our version of Murman’ s 
code and tested it. The results of a few numerical experiments quickly showed 
that it fared no better than Murman’ s second-order-accurate scheme; it failed 
to converge reliably. We believe that the reason for the failure of both of 
these second-order-accurate schemes is closely tied to the parabolic point and 
its function in stabilizing the finite-difference scheme. Consequently, before 
presenting the modifications required to make second— order— accurate schemes 
converge, we first review the purposes for using the parabolic point. 

The parabolic point in a finite-difference scheme serves the following 
three purposes: 

a. First, it excludes the possibility of expansion shocks. In other 
words, the parabolic point ensures that the fluid experiences, at most, 
a finite (and not an infinite) acceleration. This point can be seen easily 
from the transonic small-disturbance equation itself: 




+ 

^X X 


(0 ) 
y y 


0 . 


( 12 ) 
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At the sonic point, we have 


K - d) = 0 , 


( 13 ) 


which indicates that <i) =0 at the sonic point only if (1) is bounded. 

yy ^ ^xx 

b. The second purpose of the parabolic point is to ensure consistency 
in the flux conservation across the sonic line, where the switching operators 
are employed. 

c. Third (and perhaps most importantly) the parabolic point guarantees 
that the discretized system matrix can be inverted (and, if possible, constrains 
it to remain positive definite). To explain this point, we consider node 

figure 4C. At point the test of whether the node is to be treated 

as subsonic or supersonic is based on the sign of the coefficient which 

is evaluated by means of central differences. Since the point P^ is (by 
definition) supersonic, a backward difference at P^ creates an inconsistency 
in the matrix which is removed when the parabolic condition (j)^ = 0 is applied. 

All three of these conditions are satisfied in the first-order-accurate, 
backward difference scheme if a single parabolic point is introduced; however, 
the same is not true for second-order schemes. In particular, point (see 

figure 4C) again introduces an inconsistency into the matrix if a second- 
order-accurate, backward difference scheme is used there. At successive 
supersonic points downstream of point , the second-order-accurate, backward 

difference formula can be applied without difficulty. Consequently, our sug- 
gested remedy is to introduce not one, but two parabolic points at the sonic line 
when second— order schemes are used in the supersonic region. 
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Using two parabolic points, we have done numerical experiments in which 
second-order-accurate formulas were used in both the subsonic and supersonic 
regions. The results have shown that the second-order schemes do converge 
reliably (with about the same number of iterations as is required for first- 
order schemes) when the second parabolic point is added. In order to obtain 
a second-order system that is completely analogous to Murman’s first-order 
scheme, we also consider the requirements of a shock-point operator in a 
second-order scheme. 

The shock-point operator serves some analogous (though not completely 
identical) purposes in cases where the flow switches from supersonic to 
subsonic, as the parabolic point does when the flow goes from subsonic to 
supersonic. When decelerating through the sonic point (whether discontinuously , 
as across a shock, or continuously, as across a decelerating sonic line), we 
introduce the shock-point operator for the following purposes: 

a. The shock-point operator allows for a discontinuity in (j)^ . 

b. The shock-point operator ensures a consistent flux conservation in 
the presence of switching operators. 

c. The shock-point operator guarantees that the discretized system 
matrix can be inverted (and that it is positive definite, if possible). 

Note that, except for the first, these purposes are identical to those of 
the parabolic operator. Item (a) in the parabolic list ensures that (expansion) 
shocks cannot occur; item (a) in the shock-point list ensures that (compression) 
shocks can occur. 

To introduce our shock point, we start by analyzing Murman's 
shock-point operator for the first-order-accurate scheme. Using the 
nomenclature in the sketch below. 
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we have from equation (12) 


(K"l - - CKU 3 - ..^) = 0 . (14) 


Then, by adding and subtracting the quantity, 


Ku, 



Murman obtains 


(Ku^ “ i^u^) - (KU2 - + (Ku2 - hn^) ~ (Ku^ - ^u^) 


+ ~2 = 0 - 


(15) 


or 


K - h(n^ + U2) (u^ - U2) + j^K - ^(u2 + u^)j (U2 - u^) 




( 16 ) 
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Murman then linearizes equation (16) by evaluating (u^ + n^fl and 
(u 2 + u^)/2 using the most recent available values. Extensive numerical ex- 
periments have confirmed that the resulting scheme is stable and that the 
iteration converges with extremely good reliability. 

The introduction of the shock-point operator for our second-order-accurate 
scheme must also be done with care. We begin by first considering a shock-point 
operator that is analogous to Murman’ s and then consider the one suggested by 
Warming and Beam (see ref. 20). Emphasis in both of these cases will be on 
satisfying all three of the purposes outlined above. 

We use the notation in the sketch below. 


i-3, j 


-Ht 

i-2, j 


r.,. I 

I 




Y i,j+l 


i+l,j 

i,j-l 


We allow a discontinuity in (f)^ between points i-1 and i and obtain 
an elliptic shock-point operator with a derivative boundary condition for the 
velocity U 2 downstream of the shock. This boundary condition must be consistem: 
with the locally normal shock-jump relation 


U 2 = 2K - u^ 


(17) 


together, they give 


K - 



+ 

T 




U2) 






(18) 
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Numerical experiments with this scheme have indicated that it allows reliable 
convergence in the relaxation iteration and that the results converge to the 
same limit reached with Murman's first-order-accurate scheme. 

We have also obtained a second-order scheme that is formally equivalent to 
the one suggested by Warming and Beam. To obtain this scheme, we replace 
equation (17) by 

= 2K - (2u^ - u^) (19) 


and substitute equation (19) into equation (18) . Again, we obtain a stable 
scheme. Some results are given in figure 5. The effect of the shock-point 
operator is distinguishable in the fine grid calculations. 

Use of an Elliptic Solver 

In all numerical techniques for solving the transonic equation, some method 
of linearization is used to convert the nonlinear equation into a system of 
linear equations, which can then be solved by various means. The specific tech- 
nique for linearizing the equation has a considerable effect on the amount of 
computer time that will be required to obtain the solution and may even determine 
whether the iteration for the nonlinearity converges or diverges. We now attempt 
to classify some methods for linearizing the equation, with emphasis on under- 
standing how and why the Poisson technique works. 

a. Picard (linearization by freezing coefficients) 


(K - (f)^)(J)’^'''^ + = 0 

X ^xx ^yy 


( 20 ) 
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b. Rayleigh- Janzen (linearization about the equivalent incompressible 

flow) 


Trj.n+1 . ,n+l ,n ,n 

K© + © = © d) 

XX yy X XX 


( 21 ) 


c, Newton-Raphson (Let = (|)^ + 6(J> , hence J6^ = -R((})^) where J 

is the Jacobian and R(0^) is the residual. As a special case of this 
method, J(0^) is approximated by J independent of n.) 



(4)"^ + 6(J) ) + (4)^ +64) ) = 0 

XX XX yy yy 


( 22 ) 


which, neglecting second-order terms, gives 


(K - (j)^)64) + 64) - 4)^ 64) 

X XX yy XX x 


-R(4)'') , 


where the right-hand side is the residual 


R(4)^) = (K - 4)’^)4)^ + 

hr / ^x XX yy 


(23) 


Note that the Picard iteration is basically Murman’s scheme (except the linear- 
ized coefficient is updated during the iteration instead of after it) , the 
Rayleigh-Janzen scheme is the one chosen by Martin (ref. 21) and is discussed 
at length below, while the Newton-Raphson scheme is similar to Iterative 
schemes based on the unsteady small-disturbance equation. 
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In the Newton-Raphson scheme, the differential operator 


2 2 

. .nv 3,9 .n 

- ‘>’ k > - 2 -^— 2 - *XK 
dx dy 


3x * 


( 24 ) 


which is the Jacobian of the small-disturbance equation, is of mixed 
type and switches character at the same location the original equation 
does. If we represent the Jacobian as 


2 2 

1_ + + 3_ ^ 

,.2 » 2 Ax 3x * 

3x 3y 


(25) 


we note that it is very similar to the unsteady, small-disturbance equation. 


(K - 4) )(j) + 

X XX 


i - 2d) = 0 

yy xt 


(26) 


where we associate d) with d)^ - d)^ . 

xt X X 

The convergence of these iterative methods for the case of elliptic 
equations has been studied extensively; however, their application to tran- 
sonic flow represents an entirely different problem. Martin and Lomax (ref. 2) 
have suggested an iterative procedure for transonic flow, which is identical to 
the Rayleigh-Janzen technique. More specifically, they linearized the transonic 
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equation by placing the nonlinear term on the right-hand side and then solved the 
resulting Poisson equation by a direct matrix inversion Cf ast elliptic solver) . 
After each iteration, the nonlinear right-hand side was updated. Their Initial 
efforts were successful for transonic solutions that included very small super- 
sonic regions; however, after applying special stabilizing procedures, they were 
able to extend the technique to free-stream Mach numbers, which allowed much 
larger supersonic bubbles. The obvious question that this technique raises 
is: How can a mixed- type equation be solved as a series of Poisson equations? 

Jameson has reviewed this iterative approach and has shown that the tech- 
nique fails for purely supersonic flows (see ref. 3); however, with the addition 
of an additional de-symmetrizing term to the Laplacian, the iteration can be made 
to converge, even in supersonic flow. Thus, we consider the solution of the 
equation 


a6(|) + 6(j) + 6(}) = -R(d)^) 

XX Ax X yy 


(27) 


and use central differencing for ^^yy backward differencing 

for (|)^ (and R(^n)) , to obtain 


a9 (6(J)) + g(6(j)?'^^ - .) 

XX^ T/ ^1-1, j' 


+ 


3 (6<i). .) = -Ax^R? . , 

\Ay/ yy" x,j 


( 28 ) 


where 9 f = f . - , - 2f . . + f ^ . and 9 f = f . . - - 2f , . + f , _ _ 

XX 1-1, J i,J 1+1,3 yy 1,3-1 1,3 1,3+1 
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For supersonic flows with periodic boundary conditions, the Von Neumann stability 
analysis shows that stability is obtained when 3 > 2a + |Kil| , where 
Kp = K-({) ; . If 3 is chosen in this manner, we can use a fast direct solver to 
invert the left-hand side. Thus, we see that, if we attempt to solve the 
hyperbolic equation by a series of "elliptic" operators, we must add the (large) 
term 3<j>^/Ax to ensure stability, 

A heuristic analysis of this de-symmetrized operator shows that when 
3/Ax is large enough to ensure stability, the operator is no longer elliptic. 

For example, consider the following difference scheme in the x direction: 

(Ax^)L - a3^^6<(. + 6(<S<f',. - '5<t>,-_i) 

X XX 1 1*^ X 

= 0664) - 2a64)^ + 0664) + 3<S4>i - 364) , (29) 

which, after we regroup the terms, becomes 

(Ax^)L^ = (a - 6)<l)j^_j^ - (2a - g)i}.^ + a((>^^^ . (30) 

Defining the new parameter 

a' = (2a - S)/2 , (31) 

we can rewrite equation (30) as 

(Ak")L^ = (a'- - 2a'*^ +(a' + . (32) 
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which is a valid approximation for 


\ + aI [*x + (33) 

It is immediately obvious that the type of equation (33) depends upon the sign 
of cl' , which can, in turn, be controlled by the magnitude of B • Thus, if 
B > 2a , the equation becomes hyperbolic, and we see that the Von Neumann 
criterion for stability in supersonic regions is equivalent to requiring that 
the left-hand side operator be made hyperbolic. 

If we now return to equation (27) and use central differencing for (|)^ , we 
see (by again following our heuristic argument) that it is not possible to change 
the type of the equation. Similarly, the Von Neumann condition also indicates 
that the Poisson iteration for the wave equation will not converge when central 
differencing is used for the 9/3x term. Thus, we see that it is not this term 
alone, but rather the unsymmetrical differencing of it, that allows convergence. 
The introduction of the asymmetric term in backward-difference form is neces- 

sary for convergence because it removes the elliptic nature of the left-hand side 
operator and causes it to be hyperbolic. 

We have conducted a numerical experiment based on this idea by using a line- 
relaxation version of the analysis of Martin and Lomax. Instead of using the 
primitive variables u and v , we use the velocity potential ^ . The line 
relaxation should effectively introduce the precise term (J) (or 6(p ) , which 
is needed for convergence. The solution does indeed converge and the results 
are shown in figure 6. 

The special stabilization procedure referred to above, which Martin and 

Lomax used, was a u^ terra. The use of this de-sjrmmetrizing term c.llowed them 
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to extend the technique to transonic flows with large supersonic bubbles. As we 
noted above, we argue that this de-s 3 nmnetrization term is effective because 
it changes the type of the left-hand side operator in the supersonic region so 
that it is no longer elliptic. 


FINITE ELEMENTS 
Hyperbolic Methods 

Finite- element procedures can be applied directly to hyperbolic schemes; 
in particular, the Lax-Wendroff scheme applied directly over the finite-element 
formulation (finite element in space, finite difference in time). Thus, if we 
consider a Taylor series expansion in time, we have 

u(t + At) = u(t) + uAt + ^iiAt^ , (34) 


and if 


9u _ 9F (u) 
8t 3x 


we can express 


(35) 



where A is the Jacobian of F(u) , 


(36) 
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C371 


3F. 

A. . = . 

xj 3u. 


Hence, 


u(t + At) = u(t) + II At H- ^ (a II) . (38) 

Examples of using finite-element procedures to solve Lax-Wendroff schemes such 
as these have been given by Oden (ref. 22), who used a Galerkin procedure. 

If an explicit artificial viscosity term is added, we obtain an equation 
that is similar to the one studied by Wahlbin (ref. 23) . When the explicit 
artificial viscosity term is added, the transonic equation, in its small- 
disturbance form, becomes 


u = K„u + v +£u 
t Jc X y XX 


(39) 


V = u - V 
t y X 


(40) 


If we now represent the velocity in terms of a shape function Tp as 


u = ip^(x,y) U^(t) 


(41) 


V = i|;^(x,y) V^(t) 


(42) 
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and apply a Galerkin procedure, equations (39) and (40) become 

//('''i “i,t 'I'j) °ff(h Kk '(’j) '^y 



V. i(< 


,) 


dx dy + 


e ip 


l.XX 


U. 

X 


^3) 


dx dy 


(43) 


and 




dx dy = 



dx dy - 



V, Jp 


,) 


dx dy . 


(44) 


These can be written in the more compact form 


M. . 

13 


dU. 

X 

dt 


K„K. .U. + C. ,V. + ed. .U. 
£jxx jxx ij 1 


(45) 


dV, 

M. . -TT^ = C. .U. - K. .V. 
xj dt jx X jx X 


(46) 


by defining the matrices 


"If 
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ff 


( 48 ) 



ff ^"^1 

=Ji W '''j'^^ ’ 


( 49 ) 




3i|^. 

^ dx dy . 


( 50 ) 


For the special case of linear shape functions and rectangular elements, 

the functions ili. become 
1 


h 







t X y 

3 " >’l "2 





2L_\z_ 

hjh^ 


( 51 ) 


so that the matrices C, D, K, and M become 


C. . 



-2 -1 1 

-1 -2 2 

-1 -2 2 

-2 -1 1 


1 

1 

1 

2 


( 52 ) 
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D. . 



2 - 2-1 1 

-2 2 1-1 

-1 1 2-2 

1 - 1-2 2 


( 53 ) 


K. . 



-2 2 1-1 

-2 2 1-1 

- 112-2 
-1 1 2-2 


(54) 


and 


M. , 


hih2 
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4 2 12 
2 4 2 1 
12 4 2 
2 12 4 


(55) 


Implicit schemes are recommended for solving either these systems or 
ordinary differential equations in time so that stability restrictions can 
be avoided, especially in the latter case where the c{)^^ term appears. 

Kreiss and Scherer have reviewed the convergence of iterative schemes 
such as these (see ref. 24) and have shown that they are always stable for 
semi-bounded operators. Swartz and Wendroff have reached similar conclusions 
for Burger’s equation (see ref. 25). The disadvantage of using techniques of 
this type is that the mass matrix must be inverted at each step.. This 
effect can, however, be diminished by going to the lumped mass formulation. 

In passing, we mention that, instead of the Galerkin method, the least- 
squares method may be used (in space) in the same way Carasso solved the wave 
equations (ref. 26) and a coupled system of wave and heat equations (ref. 27). 


31 



We also mention an interesting work by Mote (ref 28) . Mote has considered ' 
the use of global- local finite elements, where the known characteristics of 
the global solution to the problem of interest are exploited. For example, 
the global behavior of the incompressible flow solution could be used in a 
global-local scheme for transonic flow. Note, in addition, that some global 
methods (in the classical sense of Ritz) have been applied to compressible and 
even transonic flow as early as twenty years ago by Wang (ref. 29). His 
computations did not include any shocks, although he suggested a shock-fitting 
procedure. 

Finally, we should note the important parallel between transonic flow and 
shallow-water theory. The hydraulic analogy leads to equations of motion that 
are identical to the transonic equations if the ratio of specific heat Y 
is taken to be 2.0. A comparison between finite-difference and finite-element 
techniques for shallow-water theory has been given by Weare (ref. 30). We note 
without comment that he concluded finite-difference procedures were more 
economical for this problem (see also ref. 31). 

Mixed-Type Methods 

Introductory Remarks - Some Simplified Models for the Transonic 

Small-Disturbance Equations 

According to the transonic small-disturbance approximation, the stream- 
lines are almost parallel to the x axis, and the nonlinear effects occur only 
in the x terms, as can be seen directly from the following small-disturbance 
equation: 


(K - (p )(p + <p = 0 

X ^xx yy 


( 56 ) 


32 



A one-dimensional version of this equation 


(K(j) 


X X 


0 , 


( 57 ) 


has been studied by Bauer et al. (ref. 32). By treating P as a general 
matrix and (f) as a vector, we can extend Bauer *s equation to the more 
versatile form 


(K(t>^ - + P(f. = 0 , 


(58) 


which corresponds to not one, but a system of, ordinary differential equations. 
This latter system, equation (58) , is almost identical to the transonic equa- 
tion, and it can be obtained from the transonic equation by a step that we shall 
refer to as semi-discretization. By semi-discretization, we mean that c() 
is continuous in x , but is discrete in y . As an example, we can use 
central differencing to approximate in equation (56) and so can arrive 

immediately at equation (58), where the matrix P is given by 


P 



-2 1 
1-2 1 


\ \ 

^ \ \ 

\ \ \ 
1 -2 1 


(59) 
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To complete the system, we must modify the first and last rows of P by 
applying the boundary conditions.* Similar representations can also be ob- 
tained by applying the Method of Moments or by applying a finite-element 
technique in space (in the y direction only) . 

After some appropriate discretization in y has been applied to the tran- 
sonic equation, the resulting equation (58) can be categorized for either 
subsonic flow or supersonic flow by the following two classical systems of 
equations : 

a. The case of subsonic flow is closely analogous to the two-point boun- 
dary value problem 

- (QVx + = 0 > (60) 

which is the one- dimensional, steady-state heat equation. 

b. The case of supersonic flow is closely analogous to the initial 
heat value problem 

- (M(j)) * + P({) = 0 , (61) 

which is the equation for a mass-spring system. 

A separate body of literature exists for each of these equations; tran- 
sonic flow represents a combination of both of them. In the subsonic region, 
the equation is elliptic (boundary- value in nature) , while in the supersonic 

Note that P can be diagonalized: A system of uncoupled ODE's can be con- 

structed if the eigenvalues and eigenvectors of the matrix are known. 
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region it is hyperbolic (initial- value in nature) . The nonlinearity of the 
transonic equation is essential since this nonlinearity is responsible for the 
transition from one region to the other. The nonlinearity also permits com- 
pression and expansion shock waves to occur (i.e., a discontinuous transition 
from hyperbolic to elliptic). 

The application of finite- element techniques to the two-point BVP 
corresponding to equation (60) is well established. Applications of finite- 
element procedures to initial value problems, like equation (61), have also 
been used in the field of structural dynamics; however, usually finite- 
difference methods are used for this equation. Thus, the general pattern 
is that finite-element procedures are used more often than finite-difference 
procedures for boundary-value problems, but that finite-difference techniques 
are used more frequently for initial-value problems. In the present transonic 
case, the question of which technique to use is not as obvious since the tran- 
sonic equation is of mixed-type and encompasses both types of equations. 

A Comparison of Various Discretization Techniques 

In this section we compare three distinct methods that can be used to 
discretize a partial differential equation. These three methods are finite 
differences, finite volumes, and finite elements. 

Discretization methods . - The construction of a finite-difference scheme 
is usually based on a Taylor series expansion. The stability of the resulting 
difference scheme can then be easily investigated in a linear, local sense by 
applying the Von Neumann stability analysis, as long as the grid system is 
rectangular. Heuristic stability analyses, which again rely on a Taylor series, 
can also be useful in categorizing the truncation error as either dissipative 
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or dispersive. Boundary conditions are easily included in a finite-difference 
scheme when the boundaries are rectangular, but when the boundaries become 
nonrectangular or irregular, their finite-difference representations become 
cumbersome and inaccurate. Thus, the application of boundary conditions or 
irregular boundaries represents one of the key weaknesses of finite-difference 
procedures. 

In the finite- volume technique, the principal idea is to convert the 
differential equation into its integral form before applying discretization pro- 
cedures. The advantage of this approach is that it transforms the differential 
equations, which are really mathematical expressions of certain conservation 
laws, into a form that allows the conservation of these same quantities to be 
verified and enforced easily in the discrete system. 

The third type of discretization, finite elements, can be obtained by 
applying either a variational principle (such as the one used in classical 
mechanics) or a weighted residual method. The weighted residual methods, in 
turn, include Galerkin techniques and least-squares methods. Finite-volume 
techniques can be considered as specific realizations of the method of weighted 
residuals (the method of subdomains). 

Applications of discretization techniques to elliptic problems^ . - The 
classical five-point formula for the Laplace equation represents an interesting 
example of the application of the three techniques described above (see Varga's 
discussion). For the case considered here, the same five-point formula can be 
derived by means of a Taylor series expansion (finite difference) , by conser- 
vation of flux across imaginary boundaries (finite volume) , or by means of 
finite-element techniques based on triangular elements and linear shape functions 
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as shown for the first time by Courant.* Although all these techniques lead to 
identical results for the Laplace equation on a uniform grid system, in general, 
they lead to different schemes. 

As an example of the types of differences that are generated when the 
three techniques are applied in more general circumstances, we again consider 
the Laplace operator; but this time, instead of considering the five-point 
formula, we consider the nine-point scheme in the sketch below. 



which can be expressed as the linear combination of the following two five- 
point formulae. 


Four 



Plus One 



Birkhoff and Gulati (ref. 33) have noted that Courant^s derivation of the 
five-point formula for the Laplace equation, which is based on the Ritz 
variational method, does not generalize to the Poisson or Helmholtz equation. 
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Finite- volume techniques, which are based on replacing the area integral 


by an equivalent line integral (Gausses theorem) in the following fashion. 


n 


give a result which is equivalent to the finite-difference representation 





r 


if linear shape functions are used on rectangular elements, with a weighting 
function, which is unity, inside the domain and zero outside. Note that be- 
cause of the cancellation on the interior sides of the line integration, only 
the corner points appear in the finite-volume representation. 

Finally, the use of finite-element techniques based on a variational 
principle, result in a discretization that can be expressed as 
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and that, in finite-difference terras, is equivalent to the equally weighted 
sum of the finite-difference formulae 



Birkhoff and Gulati have considered general discretization procedures for 

linear source problems, with a view towards determining optimal, few-point 

representations (see ref. 33). Their comparative study considers both five- 

point and nine-point discretizations on regular mesh, along with some 

discussion of the three-point analogs for the corresponding one-dimensional 

case. The results show that global accuracy of five-point formulae could 

never be greater than 0(h) when the grid was nonuniform, nor better than 
2 

0(h ) when the grid was uniform, regardless of how the formulae were obtained. 

This accuracy is, of course, exactly the accuracy that is obtained with standard 

five-point, finite-difference formulae. They also demonstrated that the opti- 

4 

raal accuracy for the nine-point formula was 0(h) with a uniform rectangular 

2 

mesh and that it deteriorated to 0(h) for a nonuniform mesh. This order of 
accuracy is equivalent to that achieved by the Rayleigh- Ritz method, with 
bilinear approximating functions. 

As an example of the use of higher-order elements, we consider the ap- 
plication of piecewise-continuous cubics to the discretization step. The 
dependent variables can be represented by cubics through the use of three 
distinct elements (shown in the sketch below): (1) cubics for which cnly the 

function itself is constrained to be continuous at the nodes (e.g., serendipity 
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elements) ; (2) Hermite cubics for which the function and its (two) first deri- 
vatives are continuous at the nodes; and (3) cubic splines (tensor product 
of a one-dimensional spline) whose second-derivatives are continuous at the 
nodes. 



Since cubic splines require more grid points and are restricted to rectangular 
elements, we actually compare only the Lagrange and Hermite cubics. The 
Hermite cubic formulation seems to be more accurate. Later, we give some 
numerical results for incompressible flow over a parabolic-arc airfoil to 
demonstrate this accuracy. For the sake of comparison, we also discuss a mixed 
variational principle that uses linear shape functions in (p , u(=({)^) and 
v(=(})y.) • As noted previously, such lower-order elements can be used in con- 
junction with extrapolation techniques (such as Richardson’s) to obtain 
higher-order accuracy. 

Wheeler (ref. 34) has analyzed another method of obtaining very accurate 

approximations of the flux values at particular points in the domain for the 

two-point boundary value problem. This method is based on evaluating the moment 

of a Galerkin solution of the problem, and it reduces the error in the flux 
IT 2ir 

from 0(h) to 0(h ) , Such improvements in the computation of the flux 

appear quite attractive for a velocity potential solution for which the flux is 
the principal quantity of interest. 
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Applications of discretization techniques to hyperbolic problems . - For the 
wave equation, the literature can be divided into two categories, namely, those 
algorithms that use finite-element techniques in space but finite differences in 
time, and those algorithms that use finite elements in both space and time. It 
is interesting to note that most of the work concerned with developing mathemat- 
ically rigorous proofs of the characteristics of the discretized system have 
used a combination of finite- element and finite-difference techniques, while 
most of the work using only finite elements has been more engineering oriented. 

Examples of work in the first category include the efforts of Birkhoff 
and Dougalis (ref. 35), Swartz and Wendroff (ref. 25), Swartz (ref. 36), 
Vichnevetsky and De Schutter (ref. 37), Vichnevetsky and Pfeiffer (ref. 38), 
and Goudrea and Taylor (ref. 39). On the basis of their work on the wave equa- 
tion, Birkhoff and Dougalis recommend the Numerov scheme, which is a combination 
of both finite differences and finite elements in space. The Numerov scheme 
takes advantage of the fact that the phase errors in finite-element and finite- 
difference schemes are opposite in sign, and by the use of a proper combination 
of the two phase errors, obtains a scheme with excellent dispersive properties. 

Goudrea and Taylor evaluated different numerical integration methods in 
structural dynamics, including the methods of Newmark, Wilson, and Houboult. 
Argyris et al. (ref. 40) have generalized the Newmark family of schemes and 
they were able to obtain unconditionally stable schemes by incorporating 
a few simple modifications. 

Finite-element techniques in space and time have been studied by Argyris 
(see refs. 40 and 41), Fried (ref. 42), and Zienkiewicz and Lewis (ref. 43). 

The work of Argyris and Fried is based on Hamilton's principle. In their 
work, they do not allow the function at the end point to vary C^s is tradition- 
ally done in a Hamiltonian approach) , but they do allow the magnitude of the 
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initial velocity to vary. This method is similar to the method of inverse 
shooting, which was used in our previous paper (ref. 44). In that work we also 
replaced the initial value problem by an equivalent boundary value problem. 

Zienkiewicz and Lewis have based their work on Galerkin or least-squares 
techniques. They use Hermite cubics and consider both the end position and the 
velocity as unknowns; thus, they obtain a series of weighted residual equations. 

Some algebraic examples of finite-element formulae for initial value and 
boundary value problems . - To illustrate the differences between these various 
discretization techniques, we consider a simple, one- dimensional example problem: 




XX 


- Cj) = 0 . 


(62) 


We have constructed formulae for both elliptic (boundary value) and hyperbolic 
(initial value) problems (using Hermite cubics as shape functions) for each of 
three finite-element procedures, namely, the Hamiltonian, the Galerkin, and the 
least-squares techniques. The results are summarized in figures 7 and 8, where 
we present the influence coefficients corresponding to equation (62) for the 
node between two adjacent (one-dimensional) elements. 

As an example, we describe the procedure for determining the Galerkin 
results in detail so that the diagrams in figures 7 through 9 can be more 
clearly understood. 

We define the numerical representation of (p in equation (62) over any 
element in terms of the Hermite interpolating polynomials 3^(x) , as 

$ = 3^(x)x^ , i = , (63) 
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where the parameters represent the unknowns Cthe influence coefficients) 

determined for each element. The Hermite polynomials are cubic curves that are 
defined over each element and that satisfy the following four sets of boundary 
conditions at the ends of the interval (0,h) . (Since the polynomials are 
cubic, we can specify four end conditions in each interval.) 



Thus, each polynomial satisfies one unity boundary condition and three homo- 
geneous conditions. Algebraically, these polynomials are defined as 


3j^(x) = 1 - 3(x/h)^ + 2(x/h)^ , (64) 
^^(x) = 3(x/h)^ - 2(x/h)^ , (65) 
3^(x) = ^x/h - 2(x/h)^ + (x/h)^j h , (66) 
3^(x) = l^-Cx/h)^ + (x/h)^j h . (67) 


By the usual Galerkin procedure, we choose a weighting function ip and 
require that the equations, after being multiplied by 4^ , be satisfied on the 
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average over the interval (0,h) corresponding to one element. Thus 
we require 



dx = 0 


( 68 ) 


and immediately integrate by parts to give 

= B.T. , 


( 69 ) 


where B.T, refers to boundary terms occurring at the ends of the element. 

We now choose the weighting function to be equal to the polynomial 3^ 
corresponding to the appropriate degree of freedom. Then, using the numbering 
system defined below for a single element, 


^ 1 

$ 3 

X 


-K 

2 

4 


we compute the influence coefficients, in turn, for the function and its 
derivative at the right end and, then, at the left end. 

The influence coefficients for the function (p^ at the right end of the 
element are given by 
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/oNv, + = lo"' |«i®2 + ^ie2>Xl 


+(e^^ + 62^X2 + (P3P2 + 

+(6482 + ®4^2^X4 } dx = B.T. , ( 70 ) 

where the weighting function for point 2 is 32* After performing the indi- 
cated integration, we obtain the following relation between the four values 
of Xi , 


/ 54h _ 36\ / l56h 3e\ 

\420 30/^1 \ 420 30/^2 


(l3h! _ i_\ + ^ _ j_\ 

\420 10/^3 420 10/^4 


B.T. 


(71) 


If we repeat this integration for the degree of freedom number four 
((})x at the right end) here, using 3^ p-s the x^^eighting function, we obtain 
a similar relation, and similarly, for the two degrees of freedom at the left- 
hand end. The results are summarized in figure 7 . 

Having derived the influence coefficients for each degree of freedom 
in terms of the other degrees of freedom in the same element, we now combine 
two elements and obtain the appropriate influence coefficient for the center 
point in terms of the six Ctotal) degrees of freedom in the two elements. De- 
fining two elements A and B and using the following notation. 
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A B 

X )( X 

$ 1 2 3 

4 > 4 5 6 

X 


we obtain, for the second degree of freedom (for the value of (j) at the 
center), which is the left end of element B and the right end of element A: 


54h 

420 


36\ +/i 5^+36\ /mi 

30/^1 \420 iOf^l \420 



X4 


420 30 1 ^5 \420 30 j ^2 




10/^' 


= 0 , 


( 72 ) 


which is an equation for X 2 terms of X 3 » X^> X 5 » Xg • ^ similar 

application gives a second equation for degree of freedom number five, which 
corresponds to the derivative (|)^ at the intersection between the two elements. 
The influence coefficients for ({) are summarized in figure 8A, along with 
similar results that have been obtained with the least-squares procedure. The 
results obtained from the Hamiltonian method are identical to those for the 
Galerkin procedure and have not been rewritten in figure 8A. Figure 8B 
gives the same comparison of influence coefficients, but this time for cj)^ . 

The supersonic (initial value) case can also be computed from the 
results shown in figure 7; but, when the two elements are considered together, 
the downstream element must not have any effect on the solution for the in- 
fluence coefficients corresponding to the central point. Again, using the 
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notation from the sketch on the previous page, we write the relationship for 
X 2 the supersonic case as 



36\ +/l5^+36\ 

30/^1 \420 301^2 


+ 


/ 13h^ 

3 \ . 

1 22h^ 

3 \ 

^420 

30 j ^4 ^ 

y 420 

30/ 


(73) 


These Galerkin results for the supersonic problem are summarized in 
figure 9 in a format identical to that used for the subsonic case (figure 8). 
Again, Hamiltonian and least-squares results are also shown. Note that in the 
initial value case the Hamiltonian results (where the method of inverse shooting 
has been used) differ from the Galerkin results. Also note that the downstream 
points do not contribute to the solution at the central point. 

Influence coefficients for a two-dimensional example . - Similar tables for 
a two-dimensional case (partial differential equation) are given in figures 10 
through 12. For this case, we have considered the Laplacian equation: 


(j) + 

XX 


yy 


0 


for the elliptic case, and the wave equation, 


— (b + (b = 0 

XX 


(74) 


(75) 


for the hyperbolic case. Again, these tables are based on the Hermit ian inter- 
polating polynomials. 
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Chan and Brashears (ref. 45) have applied a least-squares procedure to 
transonic flows. In this work, the upstream effects are zeroed out in the 
supersonic region. This procedure is identical with the one proposed by 
Zienkiewicz and Lewis for the wave equation (see ref. 43). 

Schemes for Transonic Flow 

In a previous paper, Hafez, Murman, and Wellford (ref. 44) derived two 
distinct finite-element schemes for the transonic small-disturbance equation. 
These two schemes used different discretization criteria for their development. 
The first scheme uses a finite-element technique in the space-like variable 
(y) and a finite-difference representation in the time-like variable (x) . The 
only difference between this finite-element scheme and Murman *s scheme is in the 
mass matrix. If a lumped mass is used, the scheme reduces identically to 
Murman* s; if a consistent (Graham) matrix is used, the schemes remain distinct. 
In the consistent-mass matrix formulation, three y levels of the x equations 
are coupled. Comparison with Murman* s results show that the two calculations 
are almost identical (differences occur only in the fourth decimal place) ; 
i.e., the lumped-mass (finite-difference) and the consistent-mass (finite 
element in space) methods give comparable results. 

The second scheme described by Hafez, Murman, and Wellford is the "inverse 
shooting" technique, which uses a finite-element discretization in both the 
space-like and the time-like variables. Unlike the first scheme, the inverse 
shooting scheme is not unconditionally stable. We can show this readily by 
writing the scheme in the form of Von Neumann and Lees: 


£ XX 


= 


C.D. 


yy 


+ h^W(l) 


C.D. 


xxyy 


K^= K-^ 


X 


(76) 
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where the subscript C.D. refers to central difference (see ref. 46). Von 
Neumann and Lees showed that this equation was unconditionally stable only 
when w ^ . The corresponding value for 03 in our scheme is O) = (Kil + l)/6. 

By analogy to the system of Von Neumann and Lees, we see that our system is 
non-dissipative and that it has dispersive properties that depend on 0) . 

Because of the conditional stability of equation (76), it is evident 
that terms proportional to ^yyx dissipative term) are 

needed. In particular, if these terms are added with the proper multipli- 
cative coefficient, the scheme becomes identical to the backward-difference 
scheme used by Murman. If linear shape functions and rectangular elements are 

used, the finite-element representation of both (1) and d) becomes iden- 

xxyy yyx 

tical to their counterparts obtained from centered finite differences. When 
higher-order shape functions are used for the finite- element representation, it 
no longer remains obvious that these two terms sufficiently guarantee stability. 
Fortunately, however, Showalter (ref. 47) has obtained a rigorous answer to 
this question. He studied two nonstandard methods for integrating the initial 
value problem 

B$ + A(f) = f (77) 


in a Hilbert space using Galerkin projection techniques. The first method, when 
applied to the equation 



f 


(where A refers to the Laplacian) , is 


(78) 


49 



( 79 ) 


- eA 1^ - A(j) = f(K,t) , 

9t 

where e is restricted to be positive. The second method, applied to the same 
equation, has the form 

- eA ^ - A(j) = f(x,t) . (80) 

In this model, these regularizations represent artificial viscosity or artificial 

2 2 

inertia. In our simple example, A is replaced by 9 /9y so that the first 

term is (b and the second term is (b . Some numerical examples, using 

yyt yytt 

artificial viscosity and artificial inertia, that have been obtained from our 
finite-element model are shown in figures 13 and 14. Again, we have used line- 
relaxation methods to obtain these solutions. 

As discussed in our previous work, the introduction of blending elements 
is crucial to the success of a transonic-flow computational technique. For 
example, if we use a centered finite-difference approximation of in the 

sonic element rather than a Galerkin approximation, the calculation sometimes 
diverges, and, even when it does converge, more iterations are required. Some 
examples of this behavior are given in figure 15. 

The second blending element used in these calculations is a f ini te- element , 
shock-point operator, whose specific form is shown in the sketch below. 
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The same results can be obtained with locally normal shock-fitting, namely 



where velocities that satisfy the normal shock polar 

relation. These results are shown in figures 13 and 14. 



Elliptic Methods 


A Mixed Variational Principle for Transonic Flow 
When developing finite-element approximation methods, we might benefit 
from a variational development of the appropriate differential equations. In 
this variational development, we introduce functionals whose Euler equations 
are equivalent to equation (1). Initially, we introduce the "primal functional" 
J ((})) where cf) is the perturbation velocity potential. This functional is 
given as follows: 




- i h • ( 81 ) 

2 2 

where = 1 - , K 2 = ^ 00 ^^ moment we disregard the boundary 

conditions. For simplicity, we consider the Dirichlet boundary conditions ((f) 
is specified on all boundaries). The first variation of the functional J((f)) 
is defined for arbitrary r\ satisfying the boundary conditions by 


6J((|)) = LIM 
e ^ 0 


J((f) 4- £T1) - J((f)) 

e 


(82) 


2 

The second variation 6 J((f>) of the functional J((fi) is defined similarly. 
By setting the first variation of the functional equal to zero and applying 
the fundamental theorem of the calculus of variations, we get the following 
result; 

Theorem 1 : The Euler equation corresponding to J ((f)) is 
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(K, - K„ (j) )(f) + (J) 0 . 

1 2 ^xx yy 


( 83 ) 


Taking the second variation of J ((})) , we get 


Theorem 2: For arbitrary ri satisfying the boundary conditions 


= /nK - S kK * i] 


(84) 


From equation (84) we determine a basic property of the variational method. For 
subsonic flow, 6 J(<f>) ^ 0 because ^ ^ strongly super- 
critical (transonic) flow, < 0 for many points in Q . Thus, 

the second variation of the functional J((|)) has no definite sign for the case 
of transonic flow. A method for iteratively approximating this problem (essen- 
tially, the variational analog of the standard procedure of shifting the 
nonlinear term to the right-hand side of the equation), however, is to find the 
critical points of a new functional J((f)(n+1)), which is defined as follows: 



2^x 


^(n+l) 


dx dy . 


(85) 


The Euler equation of this functional is 


1 ^xx ^yy 2 X xx 


( 86 ) 
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The second variation of ) is 


6 ^ 


■L 


(Kt ri^ + n^) dx dy 
V 1 'x y^ 


Since > 0 , we get 


(87) 


6^ J({j)^^^^^) ^ 0 , for all n . (88) 

2 — 

Thus, in the iterative scheme, the second variation 6 J is always positive, 

2 

while in the original problem the second variation 6 J is indefinite. This 
inconsistency serves as a mathematical (rather than physical) explanation of 
the nonconvergence of iteration equation (86) to supercritical flow solutions 
observed by Martin and Lomax (ref. 2) and Chan and Brashears (ref. 45). The 
discrepancy noted above is the motivation for the approach developed in this 
paper. 

A mixed formulation can be developed by letting (f) be the perturbation 
velocity potential and by letting u be the x component of the perturbation 
velocity. Then, a mixed function I((j),u) associated with the small- 
disturbance transonic problem is given as follows: 


I(<I>,U)=|/ (K^4>^ + (fiy)dx dy 


- ^ f K„u^(f) dx dy + ~ f K^u^ dx dy 
z z x j z 


(89) 


We define the first variation of I with respect to (}) and u by 6^1 and 
6^1 . By setting 6^1 = 0 and 6^1 = 0 , we obtain the following result. 
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Theorem 3: The Euler equations corresponding to the functional I(^,u) are 


- K.uu + (j) = 0 

l^xx 2 X yy 


(90a) 


u “ 4» =0 

X 


(90b) 


The second variations of I with respect to and u are denoted by 6^1 
2 

and 6^ I , respectively- Erom equation (89) we obtain the following. 

Theorem 4 : The second variations of the functional I(4),u) relative to the 

parameters (p and u are 




(91a) 


<5^1 = / K,<!i dx dy , 

'n 2 X 


(91b) 


where n is a variation in <j) or u . 

Clearly, the second variations behave as shown below: 




6^1 >0 <|)^ positive , 


6 I <0 d) negative . 
u X 
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The present formulation thus allows us to divide the domain of the flow H into 
regions (usually two in number) In which the sign of the second variation of the 
function with respect to and u is always known. In contrast, if we use the 
primal formulation for transonic flow, the sign of the second variation is in- 
definite so it is not possible to specify specific regions in which the second 
variation has a specific sign. Suppose we let be the part of the domain 9 . 

in which is positive and be the part of the domain in which 

is negative. Then 0 = . We let be the functional I restricted 

to . We let I^((|),u) be the functional I restricted to . Then, 

the solution to problem (1) can be characterized as the set (<|>*,u*) such that 


MIN MIN I (4),u) = T (4)*,u*) (92a) 

u (|) 

MAX MIN I (4),u) = I^(<i)*,u*) . (92b) 

u ({) 

In Theorem 3, we have verified that the correct differential equations 
[corresponding to equation (1)] result from the mixed variational principle. 

We now rewrite the functional to include the proper boundary terms. In approxi- 
mation procedures, the set of boundary conditions for the infinite domain fi 
is normally replaced by the following set on a finite domain (also called 

a ): 
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<P ~ far-field solution 


(93c) 


For simplicity, we neglect <p^ relative to 1 in equation (93a) . The resulting 
boundary condition is 



on a slit 



(94) 


This boundary condition can also be handled iteratively as in Chan and Brashears 
(ref. 45). We included this iterative application of the boundary condition 
in our analysis, but we do not discuss it here. Then, we let ds be the arc 
length along in x,y space. We now introduce a mixed functional 

I((j),u) , which includes the natural Neumann boundary condition along : 


I(c)),u) = y f (K (j)^ + (p^)dx dy 


/ K„u^(j) dx dy + ^ r K„u^ dx dy 


-/ 


3Hb 


dx 


(})(s) ds 


(95) 


Taking the first variation of I((f),u) with respect to (}) , we obtain, as Euler 
equations, equation (90a) and 


= ^ 
dx * 


( 96 ) 
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Since the functional I((|),u) is developed for the small-disturbance formula- 
tion (which assumes a very thin body) to within the order of the approximation 
used, ~ • Thus, the boundary condition [equation (94)3 is satisfied. 

When varying !(<{),u) , we obtain one additional integral associated with . 

This term occurs when we integrate the second term on the right in equation (95) 
by parts in x , and we obtain an integral over dy (the thickness of the air- 
foil) . Since in small-disturbance theory we assume the thickness of the airfoil 
is small, we can assume this term is negligible. In fact, in most finite- 
difference approximations, the boundary conditions are applied not on the body, 
but on the chord. Then, the integral in question disappears. 

A Mixed Finite-Element Model For Transonic Flow 

To develop a finite-element model for transonic flow, we divide the 

domain ^ into finite elements . Then, 

e 

E 

= [J , 

e=l ^ 

where E is the total number of elements in the domain. On each domain 

we introduce an approximation for the potential function <j) and the per- 
turbed X velocity component u of the following form: 

$^(x,y) = 'i'^(x,y)<l>^ , (97a) 

U^(x,y) = 3^(x,y)U^ , (97b) 
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where 'i'. (x,y) and 3 ^(x,y) are finite-element Interpolation (shape) functions, 
and (j) and U are values of the potential function ({) and the velocity u at 
the nodes of the element. 

The index notation (involving i) in equation ( 97 a) and ( 97 b) implies that 

a summation should be performed over the number of nodes in the element . In 

terms of (^ ,U ) , the fianctlonal I takes the following form: 
e* e 


I($ ,U ) = -- K. . ~ L. U.U.$, 

e* e'^ 2 13 1 3 2 X3k 13k 



ijk 


U.U.U, 
1 3 k 




( 98 ) 


where 


K. . 

ij 



(K, 'i'. 'F. + 'F, 'F. ) dx dy , 

1 X J 1 J 

X -^x y y 


L 





K2 3 . \ dx dy , 

X 


"•iak “ /s 


fl >^2 ®k 

e 



dx 3 


ds . 


F. 

3 
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The finite-element analog of the variational procedure [equation (92 )J is to 


require that the functional 
variations in and U. . 


3 

9l($ ,U ) 
e e 

9$. 


1 


0 , 


3 


I($ ,U ) 
e e 


We set 


have a stationary value relative to 


(99a) 


9l(<^ ,U ) 
e e 

9U. 

1 


= 0 


(99b) 


From equation (98) in conjunction with equation (99), 


K. . 

13 1 


1_ 

2 


L., .U.U. 
ikj 1 k 


+ F. 
3 


(100a) 


91 

9U 


•L. U.$ 

13k 3 


+ M. U.U, 

13k 3 k 


(100b) 


= 0 . 


Equations (100a) and (100b) represent the finite-element equilibrium equations 
for a single element . The corresponding equations for the entire domain 

are obtained with standard assembly techniques (see ref. 22). 

Dual Iterative Solution Algorithms 

The use of a combination of direct solution and gradient algorithms to 

solve the algebraic equations obtained from mixed-finite element models was 

initially proposed by Ciarlet and Glowinski in conjunction with the blharmonic 
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equation (unpublished data from P. Ciarlet and R. Glowinski, 1975). These 
methods have become known as "dual iterative methods." We adopt the same 
approach here. The variational problem [equation (92)] involves a mini- 
mization procedure in (}) and a minimization-maximization procedure (in 

different domains) in U . We solve for the (j) variable, using a direct 

2 

solution method, since 6^1 is always positive, and we solve for the U 
variable, using a gradient method, which locally accounts for the sign change 
of 6^1 . We thus introduce the following direct-gradient algorithm for the 
small-disturbance transonic flow calculation (for p = constant > 0 ) : 


K.. + F. . 

ij 1 2 ikj 1 k 3 




(n+l) ^ (n) + c . 

U U + L P gy 

J 


where 


+M.. 

dUj j ik 1 k j ik 1 k 


C = 1 


£ 0 , 
j 


C = -1 


6y > 0 , 

j 


where 
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~2 — Cti) 
2 


g2Y(n) 


-L 

jjk k 





for each element. From (100b) we see that the gradient of I relative to U. 


is zero if = 0 for all nodes in the element. In fact, we see that the 
variation of I with U for positive (p^ varies, as shown in the sketch 
below. 



To prevent convergence of the algorithm to the trivial solution U = (J)^ = 0 , 
we choose C as shown in the chart below: 


4> Positive 

X 


V 

4^ 

negative 

negative 

negative 

positive 

positive 

positive 

positive 

negative 















<p Negative 

X 



C 

negative 

positive 

1 

negative 

negative 

1 

positive 

negative 

1 

positive 

positive 

-1 


We treat the gradient solution step as 
equation in the incremental time parameter 


the integration of a time-dependent 
t 


_J 

9t 


3l(n) 


At 3U. 
3 


( 103 ) 


where At is the increment of incremental time t , or 


(n) 

^ 


Numerical experiments have shown two deficiencies in the iterative algorithm. 
First, for transonic flow, there is no convergence to shock-wave solutions. 
Second, the iterative scheme seems to oscillate about the solutions. 
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Artificial Viscosity Models 

To obtain the correct generalized solution for the transonic flow problem, 
we introduce artificial viscosity. Let be a positive constant. The gra- 
dient algorithm augmented by artificial viscosity is 


CD P^l 

U = -^ U(U-d) ) - U 
t At At X 


(104) 


We use equation (104) in regions where (f)^ is positive (C = -1) . Differen- 
tiating equation (104) , with respect to x , introducing the result in equation 
(90a), and assuming U ~ , we obtain the following equation: 


K <j) + (j) = K (f) (J) 

l^xx yy 2 X XX 


K^Ei 


XXX 


K2At 

2p 


xxt 


(105) 


If we assume that the iteration converges, the last term on the right goes to 
zero, and the resulting equation is 


1 XX 


4) = K^({) (}) 

yy 2 X XX 


K2.1 


XXX 


(106) 


where is the artificial viscosity term. This term damps the cf) 

equation. To provide convergence in the U equation, we introduce another 
viscosity term into equation (104) . For a positive constant (used 

only when > 0 ) , we get 
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^ U(U-<() ) 
t At X 




At 


U 

X 



U 

XX 


(107) 


The last term on the right introduces a viscous term into the U equation. In- 
troducing a finite-element analog for the last two terms on the right-hand 
side in equation (104) and incorporating these into equation (101b) , we obtain 
the following expression: 


u(n+l) 

j 


U. 


(n) 


+ Cp 


3Y(n) 

9U. 

J 


At 1 


At xj X 


(108) 


where 


°ij “ la^ ®ij “ • 

We developed a computer program DUALIT to solve the equations formulated 
above. This code contains all of the features of the subsonic programs, 

SUBSONl and SUBS0N2, described in the Appendix. 

Numerical Results 

We did test calculations for flow about a 6-percent parabolic-arc airfoil 

with the mesh of eight-node serendipity elements shown in figure 16. For 

-2 -3 

these calculations, = 1.4 x 10 and £2 = 0.7 x 10 . In figure 17 
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we compare the resulting distribution to the results of Murman for 

K = 1.8 (M^ = 0.875) . 

A GENERAL SHOCK-FITTING PROCEDURE FOR FINITE ELEMENTS 

A more general shock-fitting procedure, which was previously used by 
Wellford and Oden in elastodynamics (see ref. 48), can be applied to transonic 
flows. We demonstrate this application by starting from the unsteady 
transonic small-disturbance equation from Cheng and Hafez (ref. 14): 


acj)^=(K-<l))(/) +(i) 

xt X XX yy 


(109) 


We obtained this equation directly from the unsteady velocity potential equation 
by ignoring the "high frequency" term also use it to re- 

present the iterative (relaxation) procedure (except perhaps with a different 
value of a). We proceed with development of our general shock-fitting pro- 
cedure by obtaining the weak solution of this time- dependent equation, which 
will enable us to predict the shock's movement (both its displacement and 
its speed) . 

We represent the shock shape by the function 

f(x,y,t) = 0 , (110) 


which we assume can be inverted to give the x position of the shock 


X - x^(y,t) = 0 


( 111 ) 
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Using this nomenclature, we write the jump conditions admitted by the weak 
solution of the unsteady equation as 




“ ) f (f) E = 0 , 

X y y“ 


(112) 


where [[ E represents the jump across the discontinuity. With equation (111), 
we can rewrite equation (112) as 


a 



hj + 





= 0 , 


(113) 


and defining the speed of the shock as S = 8x/8t , we have 



Sa = - [[K((>^ 



(114) 


The corresponding jump condition admitted by weak solution for the irrotation— 
ality condition 




(115) 


is 


|[f(j) - f(j)E = 0 

“■ x^y y^x" 


(116) 


or 





(117) 


= 0 . 
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Combining equations (114) and (117) gives 



or, in its final form, the equation for the speed of the shock is 

^ = j<‘^- <fx> + (ir) '} • 

If the flow field approaches a steady state, the shock speed vanishes 
that in the steady state, equation (119) reduces to 

+(f)' = 0, 

which is of course the jump relation obtained from the steady-state 
equations . 

The strength of the unsteady shock II evaluated from 
Hadamard kinematical compatibility equation 


[^1 


dt 


[8F 


’9F 


■9f] 

U^t. 


3x, 


.9yi 


L. 


where the subscript S refers to the shock. For our case, F = , 


(118) 


(119) 


so 


( 120 ) 


the 


( 121 ) 


hence , 
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(122) 


(123) 


The shock jump relation can be incorporated into a finite-element procedure 
by applying finite differences in time to equation (119) in a manner analogous 
to that of Wellford and Oden: 


- x" = -B j<K - UJ + (f) ' j ” . (124) 

where 3 = At/a . (Note that we can use a similar formula for relaxation 
procedures, except that now 3 corresponds to a relaxation factor.) A similar 
finite-difference procedure can be used in conjunction with equation (122) to 
update [[(J’ U • As in most numerical procedures, an explicit technique of the 
type suggested in equation (124) imposes a maximum allowable step size because 
of stability considerations. 

As a passing remark, we note that when the transition from subsonic to 
supersonic (or from supersonic to subsonic) is smooth, the present finite- 
element algorithm is consistent since we automatically require that = 0 . 

A similar shock-fitting procedure can also be applied to the full poten- 
tial equation: 

p = -(pU) - (pv) , (125) 

I- X y 
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where, because of the irrotationality assiimption, we can express the density 
p as a function of the velocity only as follows: 


p = p(u,v) 


(126) 


Defining the velocity potential cf> as 


u = d) 

X 


V = d) 

y 


(127) 


so that 




(P 


xt 





> 


we can write the expression for the weak solution of equation (125) as 


[[pfj. + puf^ + pvfyj = 0 , 


(128) 


where 


f 

t 


: f : 


X 


f = 

y 
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Then, following the development of the small-perturbation equation, we arrive 
at an analogous expression for the shock velocity and strength. Although this 
procedure appears quite promising, it has not yet been tested in either a finite- 
difference or a finite-element formulation. 

EXTENSIONS TO THE FULL POTENTIAL EQUATION 

Although the full velocity potential equation is algebraically more com- 
plicated than the small-disturbance equation, the crucial factor that makes the 
full potential equation more difficult to solve is that the direction of the 
velocity vector in the full potential case is not known. Jameson (ref. 3) has 
developed some type-dependent finite-difference schemes and some relaxation 
techniques that work well for the full potential equation, but the corresponding 
finite-element analog is not obvious. 

Hyperbolic Methods 

One exception to this difficulty is when the equation is formulated in its 
complete unsteady form. In this case, the equation is always hyperbolic, and, 
as a result, the use of type- dependent differencing (which relies on a knowledge 
of the orientation of the local velocity vector) is not required. Consequently, 
finite-element techniques (in space) can be applied in a straightforward manner. 
This capability for bypassing type-dependent differencing and retaining central 
differences throughout the flow-field is not limited to the unsteady Euler 
formulation; other hyperbolic schemes can also be used. 

As a first example, we consider the fully hyperbolic scheme 

u = (pu) + (pv) (129a) 

t X y 
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Cl29b) 



Magnus and Yoshihara (ref. 12), using a finite-difference technique, solved this 
formulation numerically, but abandoned it because relaxation methods (velocity 
potential) were much faster. They used the Lax-Wendroff finite-difference 
scheme, with a special mesh arrangement near the airfoil leading edge. 

The weak solution consistent with equations (129a) and (129b) is 

[u]l - SpuJ - |pv]] fy = 0 • (130) 

Note that this weak solution is different from the one discussed in the proceed- 
ing section where the continuity equation 

p + (pu) + (pv) = 0 (125) 

L X y 

was used. Also, note that the flow field described by equations (129a) and (129b) 
becomes Irrotational only in the limit of a steady solution. 

As indicated above, finite-element procedures can be applied directly to 
hyperbolic schemes such as these. In particular, we can take the Lax-Wendroff ! 
scheme directly to a finite-element formulation (finite element in space, j 

finite difference, or iterative, in time) and thereby introduce the artificial 
viscosity that is necessary for the convergence of the hyperbolic formulation. 
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Elliptic Methods 


We begin by considering the full potential equation 



0 


or 


(p<f) ) + (p4> ) 

X y y 


= 0 


and identify three classical schemes that can be used to linearize the 


as follows: 

a. Picard method (linearization by freezing coefficients). 



^n+1 

n 


+ 



or 



(132a) 


(132b) 


equation 


(133a) 


(133b) 
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b. Rayleigh-Janzen method (linearization about the equivalent incompress- 
ible flow) , 




n+1 

XX 


+ (|) 


n+1 

yy 



lU 


XX 


o n n 
2u V 

2 

n 

a 





(134a) 


or 


,n+l , .n+1 

9 + 4 ) 

XX yy 


1 / n.n , n ,n 

p + p (|) 

pn y ^x^x ^y^y 


(134b) 


c. Newton-Raphson method, 


Let 


*"+1 = d," + 


(135) 


and J6(f) = -R((J)^) 


(136a) 


or 


= J(j)^ - R(({>^) 


(136b) 


where J 


is the Jacobian and 


R(4.") 


is the residual. 


The Rayleigh-Janzen scheme is the one used by Jameson (ref. 3) and it converges 
only for subsonic flows. For transonic flows, Jameson uses a two-step iterative 
procedure, in which a fast solver is followed by a line relaxation procedure. 

The first step can be calculated by finite elements. Although there is no analogy 
for the second step, we discuss some related work below. 
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Application of Optimal Control Methods to Transonic Flow 


Following Glowinski and Pironneau (ref. 49), we consider the minimi- 
zation problem. 


/ 



- u)^ + (4)y 



dx dy 


(137) 


with the constraints 


+<!)) = 
XX yy 


(K - u)u + V 
X y 


u + V 
X y 


(138) 


Wellford and Hafez (ref. 50) used a similar procedure based on a mixed 
variational principle, which allowed special treatment of the supersonic region. 
Thus, instead of constructing a functional whose Euler equation include the 
irrotationality condition and applying the continuity equation as a constraint 
[as in equations (137) and (138)] , they constructed a functional whose Euler 
equations included both the irrotationality condition and the continuity 
condition. This functional can be written in terms of 4) and u as 


1(4), u) 


= //[i 


t(k4>„ + 


2 12 

4)^) - ^ u^4> + 

^y 2 ^x 


1 3 

3^ _ 


dx dy 


(139) 


Then, again using the gradient method, we obtain the alternative scheme. 


64) = -0) 



+ 4) 


yy 



(140a) 


6u = -o}*u(u - A ) + + e„u 

^x lx 2 XX 


(140b) 
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In their work, Wellford and Hafez require 0)*u to be positive by choosing 
0)' negative when u is negative, while choosing U) to be very large, so 
that equation (140a) becomes 


Kd) + 

XX 


yy 


= uu 


X 


(141) 


Note that the term necessary to prevent the iteration from diverging. 

The artificial viscosity term will ensure proper behavior, including the 
compression shock (and excluding the expansion shock), in the supersonic region. 

These two methods can also be applied to the full potential equation, 
Glowinski and Pironneau minimize the functional 


I 


(u - (j) )^ + (v - (j) )^ dx dy 
X y 


(142) 


with the constraint 


(b + d) = w 
XX yy 


(PU)^ 


+ (Pu)y j 


+ U + V 

X y 


(143) 


and 


P 




+ v^) 


1 

y-1 


(144) 


The functional in equation (142) can also be solved by a mixed variational 
principle in terms of (f) , u , and v , by a direct extension of Wellford and 
Hafez's work. Performing the required extension leads to 
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I((j),u,v) = Jj 


(4 - u)^ + (* - v)^ 

X y 


[' 


+ p 1 (u^ - U(f) ) + (v^ - V(J) ) I - Cp^ , 

X y 


]- 


(145) 


where 


C = 1/yM' 


(146) 


The minimization of equation (145) yields the following three expressions; 


From dl/d(f) = 0 , we have 


X((j) + (}) ) = (pu) + (pv) + A u + V ; 

XX yy [_ X y J |_ x y 


(147a) 


from 9l/3u = 0 , we obtain 


2 2 

X(u-({)) = p((l) -u)+p (u(l) - u ) + p (v(b - V ) ; 

X ^x u ^x V ^y 


(147b) 


and, from 91/ 9v = 0 , we obtain 


2 2 

A(v - (f)y) = p((()^ - v) + pu(u<J)^ - u ) + pv(vc()^ - V ) 


(147c) 


where X is a free parameter. 

Additional regularization terms like those introduced for the TSDE are 
needed for the convergence of iterations based on this scheme; however, in this 
scheme, regularization terms are needed for both the u and the v equations. 
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One approach is to use, in the supersonic region, artificial viscosity terms 
analogous to those used by Jameson in his finite-difference formulation. 

For the subsonic region, the elliptic equation may be put into the Poisson 
form and solved iteratively as before. Again, convergence is guaranteed for the 
elliptic region; hence, in the subsonic region, and can be replaced 

by the simpler relations ^ ^ save computer storage space. 

SUMMARY AND CONCLUSIONS 

We have presented a review of finite-element and finite-difference 
techniques for transonic flow calculations. For the transonic small-disturbance 
equation with linearized boundary conditions, Murman's finite-difference scheme 
is preferred over finite-element procedures, especially if some of the more 
recent acceleration procedures (e.g., multi-grids) are used. Such a comparison 
with finite elements is based on rectangular elements with linear shape functions. 
If higher-order shape functions or mixed formulations with fewer elements are 
used, and if a rapid, direct inversion technique is used, the finite-element 
method becomes more attractive. Treatment of the nonlinear mixed equation must 
include the transition from one region to another, and to allow only for a 
compression shock; the technique must have some provision for excluding 
expansion shocks. 

Hopefully, finite-element methods will retain their advantages in the treat- 
ment of boundary conditions for the full potential equation. An unsteady 
(hyperbolic) approach for applying finite-element procedures to the full potential 
equation has been outlined. Without efficient methods for accelerating the 
iterative procedure, however, this method will not be economical. 
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Finally, we have looked for optimal control methods for transonic flow; 
but, without the inclusion of artificial viscosity or shock-fitting procedures, 
we do not expect the calculation of transonic flow fields with large supersonic 
regions to be successful. 

Flow Research Company 

21414-68th Avenue South 

Kent, Washington 98031 
June 14, 1978 
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APPENDIX 


FINITE-ELEMENT SOLUTIONS USING GAUSSIAN ELIMINATION 
(FRONTAL SOLVER) 

As part of this research, we developed two codes to calculate subsonic flows 
using standard finite-element techniques. Both programs solve incompressible 
and small-distrubance flows about arbitrary bodies, and both could be easily 
extended to a full potential formulation. 

The first program (SUBSONl) uses two-dimensional isoparametric elements 
of the serendipity type. 



modeled at all nodes) 

The second program (SUBS0N2) uses the tensor product of Hermite polynomials 
mapped to an irregular region with auxiliary mapping points (denoted by x) . 


80 





Y > 4 ^ » 4 ^ 


These two codes are based on the following system of routines : 



The capabilities of the routines, and, therefore, the program is indicated 
in the following descriptions: 

MESH2D - Generation of mesh of elements with curvilinear sides 
for arbitrary domains 
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PMESH 


Mesh plot 


ITER “ Iterative driver for solution of subsonic or transonic 

flow problems 

STIFF - Generation of matrix equations for each element 

ZIPP - Frontal solution (direct) routine for a symmetric matrix 

UNZIPP - Frontal solution (direct) routine for a nonsymmetric matrix 

GRAD - Computation of the gradients of the solution at the Gauss 

points 

As an example, we computed the incompressible potential flow about a 0.03 
parabolic-arc airfoil. In each case, the model contained approximately 
600 degrees of freedom. The first model had the four-node serendipity element 
with 28 elements in the x direction (x = -1.5 to x = 1.5) and 20 elements in the 
y direction (y = 0 to y = 2) • There were 12 elements on the airfoil. 

As a second example, the same case was computed with the eight-node 
serendipity element with 18 elements in the x direction (x = -1.5 to x = 1.5) 
and 10 elements in the y direction (y = 0 to y = 2) . There were 6 elements on 
the airfoil. 

The solution time for the linear algebraic equations (using the frontal 
solver) was 30.15 C.P. seconds for the four-node case and 29.44 C.P. seconds 
for the eight-node case on an IBM 370-158. This time is equivalent to approxi- 
mately 7.5 C.P. seconds on the CDC 6600 (using a factor of 4 for conversion). 
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Similar runs using line-relaxation methods to solve the same problem (as described 


in the section entitled Extensions to 
approximately 7.0 C.P. seconds for 70 
at least for linear problems, the two 
times for the moderate-sized examples 


the Full Potential Equation) have taken 
relaxation sweeps on the CDC 6600. Thus, 
solution schemes have equivalent computing 
tested here. 



REFERENCES 


1. Murman, E. M. ; and Cole, J.; AI AA J « 9 , no. 1, 1971, pp. 114-121. 

2. Martin, E. D. ; and Lomax, H. ; AIAA Paper No. 74-11, 1974. 

3. Jameson, A.: AIAA Second Computational Fluid Dynamics Conference, Hartford 

1975. 

4. Sichel, M. : Phys. Fluids 6, 1963, p. 563. 

5. Landahl, M. T. : Unsteady Transonic Flow , Pergamon Press, New York, 1961. 

6. Keyfitz, B. ; Melnik, R. ; and Grossman, B. : Grumman Report R5-525, 1976. 

7. Sirovich, L. ; and Huo, C.: AIAA J, 14 , no. 8, 1976. 

8. Landau, L, D. : Fluid Mechanics , Pergamon Press, New York, 1959. 

9. Guderley, K. : The Theory of Transonic Flows , Pergamon Press, New York, 

1962. 

10. Zierep, J. ; and Oswatitsch, R. : ZAMM 40 supp. , T143-144, 1960. 

11. Melnik, R. ; and Grossman, B. : Symposium Transsonicum II, Gottingen, 1975. 

12. Magnus, R. ; and Yoshihara, H. ; NASA CR-2186, 1973. 


13. 

Murman , E . M . 

: AIAA J. 12, 

no. 

5, 1974, pp. 

626-633. 

14. 

Cheng , H . K . ; 

and Hafez, M, 

M. : 

AIAA Paper 

73-88, 1973. 

15. 

Cheng, H. K. ; 

and Hafez, M. 

M. : 

AIAA Paper 

7551, 1975. 

16. 

Caughy, D. A. 

; and Jameson, 

A. : 

AIAA Paper 

76100, 1976. 


17. Ballhaus, W. ; and Steger, J.: NASA TM X-73082, 1975. 

18. Brandt, A.; and South, J.: Project SQUID, Monterey, 1976. 

19. Yu, N. J.; and Seebass, A. R. : Symposium Transsonicum II, Gottingen, 1975. 

20. Warming, R. F.; and Beam, R. M. : AIAA Second Computational Fluid Dynamics 

Conference , Hartford, 1975. 

21. Martin, E. : AIAA Second Computational Fluid Dynamics Conference, Hartford, 

1975. 

22. Oden, J. T. : Finite Elements of Nonlinear Continue , McGraw Hill, 1972. 

23. Wahlbin, L. B. ; Mathematical Aspects of Finite Elements , AP. 1974. 


84 



24. 


Kreiss, H. 0.; and Scherer, G,: ^J^e^tixal As pects of Finite Element in 

Partial Differential Equations , Academic Press, New York, 1974. 

25. Swartz, B. ; and Wendroff, B. : Math, of Comp. 23 , no. 105, 1969, pp. 37-49. 

26. Carasso, A.: Math, of Comp. 29 , no. 130, 1974, pp. 447-463. 

27. Carasso, A.: Math, of Comp. 28 , no. 127, 1974, pp. 757-767. 

28. Mote, C. D.: Int. J. Num. Methods Eng. 3, no, 4, 1971, pp. 565-574. 

29. Wang, C. T. : J. Aero, Sci. 17 , 1950, p. 343. 

30. Weare, T. J. : Comp. Meth. in Appl. Mech. and Engr. 7, 1976, pp. 351-357. 

31. Kawahara, M. : Second International Symposium on Finite Element Methods in 

Flow Problems, 1976. 

32. Bauer, F.; Korn, D. ; Jameson, A.; and Garabedian, P. : Supercritical Wing 

Sections , Springer Verlag, New York, 1975. 

33. Birkhoff, G.; and Gulati, S.: SIAM J. Num. Analysis 11 , no. 4, 1974, 

34. Wheeler, M. F. : SIAM J. Num. Analysis 11 , no. 4, 1974. 

35. Birkhoff, G. ; and Dougalis, V. A. (edited by Vichnevetsky , R.): Advances in 

Computer Methods For Partial Differential Equations , AICA, 1975. 

36. Swartz, B. (edited by Vichnevetsky, R.): Advances in Computer Methods For 

Partial Differential Equations , AICA, 1975. 

37. Vichnevetsky, R. ; and De Schutter, F. : Advances in Computer Methods For 

Partial Differential Equations , AICA, 1975, p. 46. 

38. Vichnevetsky, R. ; and Pfeiffer, B. : Advances in Computer Methods For 

Partial Differential Equations , AICA, 1975, p. 53. 

39. Goudrea, G. L. ; and Taylor, R. L. : Comp. Meth. in Appl. Mech. and Engr. 2, 

1972, pp. 64-97. 

40. Argyris, J. H. ; Dunne, P. C.; and Angelopoulos , T. : Comp. Meth. in Appl. 

Mech. and Engr. 1973, pp. 203-250. 

41. Argyris, J. H. ; and Scharpf, D. W. : T he Aero. J. of the Roy. Aero. Soc. , 

1969, pp. 1041-1044. 

42. Fried, J. : AIAA J. 7, 1170, 1969. 

43. Zienkiewicz, 0. C. ; and Lewis, R. W, : J. of Earthquake Engr. , 1972, 

pp, 407-408. 


85 


44. Hafez, M. M, ; Murman, E, M. ; and Wellford, L, C, : Second International 

Symposium on Finite Element Methods in Flow Problems, 1976. 

45. Chan, S. T, K. ; and Brashears, M. R. : AFFDL-TR-74-11 Wright Patterson Air 

Force Base, Ohio, 1974. 

46. Von Neumann, J. ; and Lees, M. : J. Soc. Ind. Appl. Math. 10 , 1962, p. 610. 

47. Showalter, R. W. : SIAM J. Math. Analysis 1, no. 1, 1970. 

48. Wellford, L. C.; and Oden, J. T. : J. of Comp. Physics 19 , no. 2, 1975, 

p. 179. 

49. Glowinski, R. ; and Pironneau, 0.: Second International Symposium on Finite 

Element Methods in Flow Problems, 1976. 

50. Wellford, L. C.; and Hafez, M. M. : Second International Symposium on 

Finite Element Methods in Flow Problems, 1976. 


86 






[•>] 




Figure 2. Murman’s 


E = Elliptic Point 
P = Parabolic Point 



Conservative Finite-Difference Scheme 


























Figure 4C. Second-Order-Accura 
Finite-Difference S 



Fully Conservative Scheme for 
;ion of Transonic Flow 




Case 

Shock- 

Point 

Operator 

Far 

Field 

Parabolic 

Point 

Grid 

Size 

A4 

Second-Order 

No 

Two 

Std 

B4 

Second-Order 

No 

Two 


Dx/2 
By/ 2 

C4 


Locally Normal 
First-Order 

No 

Two 

Std 

D4 


Locally Normal 
First-Order 

No 

Two 


Dx/2 
By/ 2 

E4 

First-Order 

No 

Two 


Dx/ 2 
By/ 2 

F4 

First-Order 

No 

Two 

Std 

G4 

No 

No 

Two 

Std 

H4 

No 

No 

Two 

Std 

14 

No 

No 

Murraan 

& Cole (1971) 

Std 

J4 

Warming & 
Beam (1975) 

No 

Warming 
& Beam (1975) 

Std 


Figure 5A. Second-Order Finite-Difference Solution of Transonic 
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Subsonic Case (BVP) ; Coefficients for (f) 
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Figure 8A. Comparison of Influence Coefficients For the (j) Term (Degree 
of Freedom Number Two) For the Galerkin, Hamiltonian, and 
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Subsonic Case (BVP) ; Coefficients for (j) 
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Figure 8 B. Comparison of Influence Coefficients For the (j) Term (Degree 
of Freedom Number Five) For the Galerkin, Hamiltonian, and 
Least-Squares Approaches (Nomenclature Is Defined in the Text) 
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Supersonic Case (IVP); Coefficients for 
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Figure 12A. Influence Coefficients For Galerkin Formulation — 
Supersonic Case For (j) Term 




































Galertcin Influence Coefficients For Term — Supersonic 




7 

1 

A I 

6 

■ 

R 

1 

in 


1 7 

11 

1 A 

15 

16 

1 1 

18 

Conlrlbu- 
1 ion 
From 
F. lenient 
A 

— 

a^li 

ficT 

nil 

>|20 

1 -etjh 
60 

III) 

210 


1 '^j 

‘ i(T 

13h 
A 20 

U) 

+ Lil! 

710 



1 

1 

180 

h' 

lOS 

. ? 
h 

."i 

ini 


no 

- , 2 
2'ij b 

AS 

h' 

lOl 




j 

Conir 1- 
hul inn 
From 
Element 
R 







No Con 

ribiil in 

1 From E 

lement II 








■ 

Coot r 1 1)11- 
l Inn 
From 
Element 
C 




™3*’ 

10 

+n'i 

A20 

-0.31, 

30 

+ llh 
210 


a^h 
^ ' 

I3h 
A 20 

60 

lib 

’710 





-^1'' 

90 

b' 

1 


, 2 

-rt^b 

"Tb?' 

h' 

105 

, 2 
rt^h 

' A 5 

4 

■ 

Coiil r 1- 

hut Inn 
From 
Element' 
1) 







No (;onl 

r Ibiit .loi 

From E 

Rinent D 










Note: Influence Coefficient For Each Point Is the Sum of All Terms In That Column 


'<,13 


1,10 7,11 1,17 


8.17 

c 

9,18 

D 

5, lA 

A 

6,15 

B 


Figure 12B. 


Influence Coefficients For the Galerkin Formulation — 
Supersonic Case For Term 


o 





























Velocity (U) 



Half-Chord 



Velocity (U) 


s 


.3 


• B2 Galerkln With Artificial Viscosity 
X El Finite-Difference 



Figure 13B. Comparison of Finite-Element and 



TE 

J I 1 I I I I 

.6 0.7 0.8 0.9 1.0 1.1 

;ance (x) 


e-Dlfference Results For Transonic Flow 








113 



Note: All Cases Gave Results Identical Within Three Significant Digits 


Figure 14B. Identification of Finite-Element Solutions Using Shock-Point Operator 
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Figure 15. Effect of Different Parabolic Blending Elements on Convergence 
in a Finite-Element Scheme 
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Figure 16. Grid System Used 
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